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Abstract 



This paper is devoted to the characterization of an extended family of CARMA (continuous-time autoregressive 
moving average) processes that are solutions of stochastic differential equations driven by white Levy noise. These are 
completely specified by: (1) a set of poles and zeros that fixes their correlation structure, and (2) a canonical infinitely- 
divisible probability distribution that controls their degree of sparsity (with the Gaussian model corresponding to the 
least sparse scenario). The generalized CARMA processes are either stationary or non-stationary, depending on the 
location of the poles in the complex plane. The most basic non-stationary representatives (with a single pole at 
the origin) are the Levy processes, which are the non-Gaussian counterparts of Brownian motion. We focus on the 
general analog-to-discrete conversion problem and introduce a novel spline-based formalism that greatly simplifies the 
derivation of the correlation properties and joint probability distributions of the discrete versions of these processes. We 
also rely on the concept of generalized increment process, which suppresses all long range dependencies, to specify an 
equivalent discrete-domain innovation model. A crucial ingredient is the existence of a minimally-supported function 
associated with the whitening operator L; this B-spline, which is fundamental to our formulation, appears in most of 
our formulas, both at the level of the correlation and the characteristic function. We make use of these discrete-domain 
results to numerically generate illustrative examples of sparse signals that are consistent with the continuous-domain 
model. 

I. Introduction 

In our companion paper, we have set the foundations of a general innovation framework that leads to the 
specification of a broad class of continuous-time stochastic processes (TJ- The powerful aspect of the formulation 
is that it unifies the classical theories of stationary Gaussian processes j2}, on the one hand, and Levy processes on 
the other |3J, the idea being that these processes can all be generated by applying a proper integral operator (L _1 ) 
to some admissible white noise process. We have also shown that switching to a non-Gaussian excitation (within 
the class of admissible solutions) necessarily induces a sparse behavior. This suggests that this type of modeling 
is highly relevant for modern signal processing, which is presently very much focused on sparse signal recovery. 
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Fig. 1. Innovation model of a generalized stochastic process. The process is generated by application of the inverse operator L _1 to a 
continuous-domain white noise process w. The generation mechanism is general in the sense that it extends to the complete family of (non- 
Gaussian) noises w = W that formally correspond to the weak derivative of some classical Levy process W(t). Gaussian processes are 
recovered by taking W{t) to be the Wiener process (a.k.a. Brownian motion). The output process s(t) is stationary iff. L —1 is shift-invariant. 

TABLE I 

Typology of continuous-time stochastic processes 





Gaussian 


Sparse 


Stationary 


classical ARMA theory 


Non-Gaussian 






CARMA processes 


Non-stationary 


Brownian motion 


Levy processes 




and present extensions 


and present extensions 



While the proposed generation mechanism is remarkably simple conceptually, it is not quite as straightforward to 
formulate rigorously because the underlying innovations (admissible white noise excitations = Levy noise) can only 
be properly defined in the sense of distributions |4), |5). Statisticians usually work around the difficulty by defining 
processes through stochastic integrals (Ito calculus) which avoids the explicit reference to white noise |6), (7); the 
downside of this widely-used framework is that it partly hides the system-theoretic aspects. 

The innovation model described in Fig. 1 is attractive to engineers because it establishes a direct link between 
stochastic processes and linear system theory. It also suggests that it is possible to transpose some standard 
deterministic techniques (e.g., determination of impulse responses, filtering, sampling of signals, cardinal spline 
interpolation) to the stochastic setting, which is mostly what this work is about. In other words, once one has gone 
through the effort of properly defining and understanding the notion of a continuous-domain white Levy noise, the 
remaining characterization problem can be addressed by relying on the powerful (deterministic) tools of functional 
and harmonic analysis. The non-trivial aspect is that one needs to resolve some instabilities (in the form of singular 
integrals), both at the system level to allow for non-stationary processes, and at the stochastic level because the 
most interesting sparsity patterns are associated with unbounded Levy measures (cf. JT] Section II. D]). 

In the present paper, we investigate the discrete-time implications of the theory for the extended class of 
continuous-time processes which are ruled by ordinary differential equations (cf. the typology of processes shown 
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in Table 1). The stationary Gaussian members of the family are well studied and play a central role in traditional 
system modeling, signal processing and control theory (2j, J8|, (9j. There is also a well-known discrete connection in 
the sense that the sampled version of a Gaussian ARMA process is itself a discrete ARMA process with the discrete 
and continuous-domain poles being related by the exponential map: {z n — e a "}^ =1 j8j, flo) , JIT) . Less obvious 
is the determination of the MA component of the discrete model which is jointly dependent upon the continuous- 
domain poles and zeros |12|. Another classical instance is provided by the Levy processes, including Brownian 
motion, which are commonly used in financial mathematics |13|, fl4) . Levy process are especially interesting in 
that context because of their ability to replicate jumps in price assets fl4) , fT5) . They are not as popular in signal 
processing circles, probably due to the fact that they are non-stationary; yet, it has been pointed out recently that 
they are actually very relevant because they are the processes for which some of present sparsity-based algorithms 
(e.g., TV-denoising) are statistically optimal [16|. The final important subclass is made up of the so-called CARMA 



processes — the non-Gaussian extension of the classical ARMA processes 1 17 1. Special instances of such stationary 



processes have been applied to financial modeling [18| and, to a lesser extent, signal processing p9)-pT|. 

In the sequel, we present a systematic characterization of the sampled versions of these processes. The primary 
contributions along the way are : 

• An addition to the non-stationary branch of the CARMA family via the introduction of generalized boundary 
conditions and "regularized" inverse operators for the solution of unstable stochastic differential equations 
(SDE). 

• The explicit 2nd-order description of CARMA processes and their non-stationary variants. In particular, we de- 
rive an exponential spline interpolation formula that relates the continuous and discrete-domain autocorrelations 
(resp., power spectra) of these processes. 

• The uncovering of the fundamental role of the exponential B-splines in the statistical characterization of the 
CARMA processes. Not only do such B-splines correspond to the autocorrelation function of the generalized 
increment processes, but they do allow for a remarkably concise description of the joint characteristic functions 
of the discrete versions of these processes. 

• The derivation of the discrete counterpart (finite difference equation) of the continuous-domain innovation 
model. The proposed formulation also extends to the non-Gaussian and/or non-stationary variants of these 
processes. 

The paper is organized as follows. In Section II, we briefly review the general innovation model which specifies 
the broadest possible class of continuous -time linear stochastic processes. We also recall the inverse-operator 
method of solution which results in a complete characterization of the generalized CARMA processes (Tj. In 
Section III, we investigate the discrete-domain aspects of the theory by considering the sampled versions of these 
processes. In particular, we establish exponential spline-based interpolation formulas that connect the discrete and 
continuous-domain correlations of the CARMA processes, including results for their non-stationary variants (e.g., 
Levy processes and their extensions). We explicitly determine the ifth-order characteristic function of the samples 
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of the corresponding generalized increment processes, which are maximally decoupled. This naturally leads to the 
specification of some equivalent discrete-domain ARMA-type innovation model. In Section IV, we use those results 
in conjunction with exponential spline calculus to develop numerical algorithms for the generation of CARMA 
processes with a special attention to the non-Gaussian, non-stationary scenarios. We conclude the paper with the 
presentation of illustrative examples of sparse processes in Section V. 

II. Review of continuous-time results 

We start with a brief review and discussion of the key results of our theory of generalized stochastic processes 
[1|. We also provide a summary of the notations in Table II. 

A. Generalized innovation models 

The continuous-time stochastic processes s(t) under consideration satisfy the general innovation model in Fig. 
1. They correspond to the solution of the (linear) operator equation 

Ls = w, (1) 

where the driving term w is a continuous-domain white noise process. The model has the ability to generate Gaussian 
processes, as well as a broad variety of sparse processes, depending upon the type of excitation noise. The delicate 
aspect is that the underlying noise processes w do not admit a standard (pointwise) interpretation as functions of t 
because they are highly singular. They can only be properly specified as distributions (a.k.a. generalized functions). 
Thus, the correct interpretation of ([T} is in the "weak" sense of distributions: 

(<p,Ls) = (<p,w), for all ip e S 

where the equality must hold true for any smooth and rapidly-decreasing test function tp in Schwartz's class S. The 
guiding principle is that, for any given ip, the scalar product (or linear functional) (ip, w) is a well-defined scalar 
random variable no matter how rough the actual noise process w is. 

As for the class of admissibl^] input noises, we have pointed out that each brand is uniquely characterized by a 
canonical infinitely divisible distribution pid(x) (or, equivalently, a Levy measure) which specifies the PDF of its 
"pixelated" observation (through a rectangular window) x — {w, rect(- — to)) which is i.i.d. and independent upon 
to (stationarity). 

The above innovation model is exploitable only if the whitening operator L has an inverse that is well-defined 
over an appropriate subset of S' (the space of tempered distributions). The equation is then solved formally as 

s = ~L~ l w & VtpeS, (<p,s) — (ip,L~ 1 w) = (L _1 *(/3, w) (2) 

'A stochastic process is called white noise iff. it is stationary and independent at all points. In our framework, this is equivalent to requiring 
that the random observation variables x\ = ((p\,w) and X2 = (ip2,w) are: 1) identically-distributed whenever <p2(t) = </?i(* — *o) f° r an Y 
to S K (translated observations), and 2) independent whenever ipi X ip% = (observation windows with disjoint support). 



September 1, 2011 



DRAFT 



5 



TABLE II 
Summary of notations 



Symbols 


Description 


Defining formula 


White noise parameters: 






V 


Levy measure 


/ min(l,a 2 ) V(da) < oo 


v(a) 


Levy density 


v(a) > and v(a) da = V(da) 




Levy function 


Levy-Khinchine formula 


Po(o) 


Poisson amplitude distribution 


p a (a) > and f R p a (a) da = 1 



Stochastic differential equations: 
L 

PL 

Id 

D 

P a 
N 
no 
a 

Pjv(C) = -P«(C) 

Pa =Pjv(D) 



D 



_d_ 
dt 



whitening operator 
Green function 
identity operator 
derivative operator 
first-order operator 
order of differential system (number of poles) 
order of unstability (number of imaginary poles) < no < N 
vector of poles 
characteristic polynomial 
iVth-order differential operator 
rational transfer function 



Ls = w. white noise 
L{pl} = 5: Dirac impulse 



P a = D - aid 



a = (cti,..., ajv) 

JV 

Pn(0 = C N + ■ ■ ■ + aiC + a = n (C - a n ) 

n=l 

P a = + • • • + aiD + a Id = P ai • • • P QJV 



PnU") 



Exponential B-splines: 



D a (z) 
/8a (t) 

/8l(*) 



first-order difference operator 
iVth-order difference operator 

localization filter 
exponential B- spline 
generalized B-spline 



A«/(t) = /(t) - e°f(t - 1) 

JV 

A a /(t) = A Q1 • • • A ajv /(t) = ^ d a [fc]/(t - fc) 



JV 



n=0 
JV 



n— 1 n=0 

joi - a„ / 27r 



/3l(0 = A aPL (t) = / 



D a {en, ut du> 
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where we are using a standard duality argument to move the action of the inverse operator (via its adjoint L -1 *) 
onto the test function tp. We have shown (T| Theorem 3] that a sufficient condition for this method of solution to 
yield a well-defined stochastic process s is 

for some constant C and p > 1, which puts some mathematical constraints on the class of admissible operators 
and excitation noises. The implicit requirement is that the excitation noise is p-admissible, which is a property that 
can be checked on the characteristic function of pid (x) (cf. Definition [T] Section II-C I. 



B. Nth-order stochastic differential equations 

We have demonstrated that the above operator method could be deployed for finding the solutions of the complete 
class of linear stochastic differential equations of the form 

JV M 

Y, ««D™S = h mV m W (4) 

n— 1 m— 1 

with N > M, where a n and b m are arbitrary complex coefficients with the normalization constraint a at = 1, 
irrespective of any stability considerations. The driving noise w, which constitutes the input of the system, is 
assumed to be white by default. The output s(t) is our generalized stochastic process whose sample values are 
generally well-defined due to the smoothing effect of the inverse operator L -1 . The characteristic polynomial of 
the underlying iVth-order system with Laplace variable £ e C is 

N 

Pn(C) =C N + Giv-iC"" 1 + • • • + a = J] (C - «*) = P a (0, (5) 

n=l 

and is also specifiable in term of its (complex) roots; these are collected in the vector of poles a — (a±, . . . , ajv) 
with the understanding that the notations Pn(C) an d Pa(C) are equivalent. 

The linear system specified by |4]) is causal-stable iff. all its poles are in the left complex half-plane. Under this 
classical assumption, its impulse response pi.it) — L _1 {<5}(t) is exponentially decaying. It is obtained by taking 
the inverse Fourier transform of the rational transfer function 

n (, A — QmCH _ h IIm=lC7"-7n) _ 1 ( ~ 

^ M -MM- b '< u u^-°,)-i^y " 

where Qm(C) = ^>mC AI + ^M-iC M_1 + ' ' ' + &iC + &o is a polynomial of degree M < N. The roots of Qu(0 stre 
the so-called zeros: 7 = (71, ... , 7a/)- The solution (output of the system) is then given by s(i) = * w)(t) and 
is stationary by construction (because of the shift-invariant filtering). When the excitation is Gaussian, one obtains 
the conventional continuous -time ARMA processes, but one can also generate a large variety of sparse counterparts 
of these processes by switching to appropriate types of non-Gaussian Levy noises. 

Remarkably, the proposed framework can also handle the unstable scenarios, the general rule being that each 
pole located on the imaginary axis induces one degree of non-stationarity. Our extended formulation requires a 
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special ordering of the poles where the n purely-imaginary roots (if present) are coming last. This gets translated 
in the following representation of the characteristic polynomial |5]): 

/ N—no \ / n \ 

p a {i<») = n cjw - ««) n o*w - mo (?) 

V n=l / Vrre=l / 

with ajv-no+m = J w m and ui m 6 R. It allows us to write the factorized version of the differential equation Q: 

(P Q1 • • • P QN _„ )(P^ 1 • • • P jUno ){s} = Qa/(D){u,} (8) 

where P Qri = (D — a„Id) is the operator counterpart of the Fourier multiplier (jui — a n ) and Qm(D) = 
IZm=i ^mD m - Each component P Qn with Re(a n ) ^ has a stable linear shift-invariant (LSI) inverse P^*, which 
is either causal or anti-causal depending of the polarity of a n . The only delicate step in solving (|8]l is the inversion 
of the second operator factor on the left which is ill-posed. Our contribution has been to propose a stable inversion 
mechanism that makes use of some "regularized" left inverse of Pj U(r The canonical solution is 

r „ /e JW< - e 3LUot \ do; 

w/w = / m ■( -s ^- (9) 



which, in accordance with ( 32 1, forces the output signal to vanish at t = 0. This ultimately yields the global inverse 
operator 

L" 1 = U ,« ■ - W P^„ n ■ --P^QmCD), (10) 

s ^ ' v v ' 

shift-variant LSI part 

to be substituted in the latter imposes the no boundary conditions on the output 

s(0) = 

(D-^ no Id){ S }(0) =0 ^ 

(D- JW2 Id)...(D-^ no Id){ S }(0) = 0. 
We have shown that this method of solution yields a generalized CARMA process s = L~ lr w that is mathematically 
well-defined. Such processes will exhibit a uq degree of non-stationarity due to the lack of shift-invariance of the 
elementary inverse operators l Um ,s- While the above inversion method is uniquely tied to the boundary conditions 



( 1 1 1, it is not the only possible approach. In the appendix, we show that one can impose other boundary conditions (in 
the form of uq generalized linear constraints: (s, cp m ) =0, m = 1, . . . , no), while retaining the required functional 
properties of the corresponding inverse operators I WmjVm and their adjoint. 

The simplest example of unstable scenario is Ds = w, which corresponds to a single pole at the origin: a\ — 
ju>i = and N = no = 1. Interestingly, the solution s(t) = Io t $w(t), which enforces the boundary condition 
s(0) = 0, perfectly maps into the Levy processes, although these are usually described quite differently |3j, p2) . 
The interest here is that we are constructing the Levy processes as the (unstable) limit of the non-Gaussian AR(1) 
family. We will see that this novel point of view facilitates the transposition of standard signal processing techniques 
to the non-stationary/non-Gaussian Levy setting, including the higher-order extensions of such processes. 
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C. Characteristic functional 

Under the assumption that the whitening operator L admits an inverse that meets the stability condition fi), 
we have shown that the generalized stochastic process s(t) satisfying the innovation model ([!]) is completely and 
uniquely characterized by its characteristic functional (cf. (T| Theorem 3]): 

= cxp ( / /(L~ l4 V(*)) dt\ (12) 



under the constraint that the so-called Levy function f(ui) is p-admissible for the same p as in (SV 

Definition 1: f(ui) is a p-admissible Levy function for some p > iff. (?) it admits a Levy-Khinchine represen- 
tation (cf. (T] Eq. (3)]) with some Levy density V{a), and (//) |/(cj)| + \u\ < C\cu\ p . 

In light of the argumentation that has been provided in pi, a more operational criterion that is equivalent to 
requirement (i) (by Bochner's theorem) is that the Fourier-domain integral 

«d(*)= f efM-***-^ 

JR 27T 

results in a well-defined (non-negative and normalized) probability density function. The latter is precisely the 
canonical infinitely divisible distribution of the noise (observation through a rectangular window) to which we had 
referred to in Section Hl-AI 

The powerful aspect of the formulation is that the functional Z s (ip) : S — > K, which is the conceptual equivalent 
of an infinite-dimensional characteristic function, condenses all the statistical information about the process. The 
underlying principle is that the inverse operator L _1 (generalized shaping filter) specifies the covariance structure 
(or generalized spectrum) of the process s, while the Levy function f(ui) fully embodies the statistical properties 
of the input noise w. 

In particular, the Levy function /Gauss(w; <To) = — |w| 2 (Tq/2, which is p-admissible with p = 2, characterizes 
the complete family of zero-mean Gaussian processes. It is associated with the canonical Gaussian PDF: Pid(x) = 
_L— g-^ /( 2 oo) and corresponds to a white Gaussian input noise excitation with variance ctq. The Gaussian 
processes are part of the larger symmetric-a-stable (SaS) family that is specified by f a {uj\a) — —^-\uj\ a with 
exponent a £ (0, 2] and scale parameter a > 0. One can readily verify that the latter is a-admissible. The SaS 
processes for a < 2 have the intriguing property that their second-order moments are unbounded (fat tail behavior) 
(23). 

Another important class is provided by the generalized Poisson processes, which are intrinsically sparse by 



construction [ 16 . These are specified by the Poisson-Levy function /p i SS on( w ; A,p a ) = A J (e ja " — l) p a {a) da, 
which is p-admissible for p = 1, under the assumption that <?{|a|} = J R \a\p a (a) da < +oo. The underlying 
excitation noise is made up of randomly scattered Dirac impulses with Poisson parameter A > (average number 
of singularities per unit time) and amplitude PDF p a (a). 
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D. Unified treatment of stationary vs. non- stationary processes 

The simplest scenario is when L _1 is shift-invariant with exponentially-decaying impulse response so that its 
adjoint provides a continuous mapping from S into itself. In such a case, s(t) — (pi, * w)(t) is obtained from w 
by convolving it with ph{t) = L~ 1 {S}(t) £ Li(R) (impulse response of the system = Green function of L). This 
allows us to rewrite the characteristic functional as 



Zs(<p) = exp(7 f((p£*<p)(t)) dt) 



where p^(t) = p-L(-t) (the impulse response of L -1 *) is the time-reversed version of pl- The main point is that 
Pl * <P G S so that the determination of the statistics of s(t) reduces to a white noise analysis with a modified 
test function: Z s (ip) = Z w (pY j * cp). The direct implication is that the output process s(t) is stationary and that its 
autocorrelation function (resp., power spectrum) is proportionarl to (pY *Pi)(t) (resp., - 1 — ). The argument 
remains valid for BIBO-stable systems in general because the boundedness condition Q is satisfied for all p > 1 
as long as p^ E £i(R) (by Young's inequality). 

While the formalism also lends itself to the handling of the nonstationary scenarios, the calculations become 
more involved when L _1 looses the convenient shift-invariance property. Fortunately, in the case of the TVth-order 
differential system, we have shown that there is an easy way out that leads to a unified characterization with maximal 
decoupling, irrespective of any stability consideration. Given the poles of the system a, one defines the following 
spline-related entities: 

1) the (causal) exponential B-spline with parameter a. (poles) 

W*)=/e--fl(if^)f <B) 

Jw „ = i \ J w - «» / 2?r 

2) the TVth-order (causal) finite difference (or localization) operator 

A a f(t) = A ai ---A aN f(t) (14) 

which is obtained by cascading the first-order operators A an f(t) = f(t) — e Q " f(t — 1) with n = 1, . . . , N. 
The important point for our argumentation is that the exponential B-spline f3 a (t) is causal, compactly-supported of 
size N with bounded derivatives up to order (N — 1) p4) , p5) . This, in turn, implies that the generalized B-spline 

Mt) = Q M (D)p a (t) (15) 

has the same support as (3 a (t) and is bounded as well (because Qm(D) is a differential operator of order M < N). 
The fundamental property of this latter function is that it can also be expressed as 

/3 h (t) = A aPL (t) = AaL-H^K*) 

2 When making statements about coiTelations, we are implicitly assuming that the variance of the input noise Oq is finite, which is true in 
most cases. For the processes with unbounded variance (e.g. SaS with a < 2), the linear coupling induced by pl remains the same, but we 
cannot relate it anymore to the autocorrelation function nor the power spectrum which are no longer well-defined. 
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where ph(t) is a Green function of L such that L{p^}(t) = 6(t), irrespective of the location of the poles in the 
complex plane. Putting these results together, we find that 

Z AaS (<p) = expUj((^*cp)(t))dt\ (16) 



which is a form that is much more favorable than (12i because the generalized B-spline /3l (which is compactly 
supported) is much better localized then p^ (not to mention the unstable scenario where the Green function exhibits 
polynomial growth). This last formula suggests that, instead of s(t), one should work with the generalized increment 
process A a s(t), which is necessarily stationary and much less correlated than the original process from which it 
is derived (cf. 1 1 , Properties 3 and 5]). 

III. Connection with discrete-time stochastic processes 

We will now show that there is an elegant connection between the continuous-time and discrete-time formula- 
tions of stochastic processes which is analogous to the connection that can be drawn between the corresponding 
deterministic linear system theories (25), p6) . The story in a nutshell is as follows: continuous-time processes are 
ruled by differential equations, while their discrete counterparts are solutions of difference equations. The equations 
and correlation structures are linked functionally through some generalized compactly-supported B-splines. The use 
of these B-splines also greatly facilitates the transposition of the methods of solution from one domain to the other. 

A. Discrete-domain notations 

Discrete processes and sequences are indexed using square brackets (e.g, s[k), h[k]) to differentiate them from 
their continuous counterparts (e.g., s(t) and h(t)). A sequence h[k] of slow growth (i.e., h[k] does not grow faster 
at infinity than a polynomial of k) is characterized by its z-transform H(e? u ) = X^fcez h[k]z~ k , which yields the 
discrete-time Fourier transform for z — . If h[k] = h(t)\ t=k is the sampled version of the continuous function 
h(t) with sufficient decay, then one can relate their discrete and continuous -time Fourier transforms using Poisson's 
summation formula: H(e? u ) = J^nez M w + ^ nn )- 

The localization operator A a (cf. ( fT~4| > and [1 Property 3]) is transferable to the discrete domain; its discrete 
impulse response, denoted by d a [k], is the inverse Fourier transform of D^e^) = rL=i(l — e"™ - ^), which 
coincides with the frequency response of the continuous-domain operator. The corresponding discrete notation is 
A a s[fc] = (d a * s) [k] = X^nez ^a[ra]s[fc — n], where the use of the square brackets indicates that the convolution 
operation is discrete. 

B. Sampled processes 

Here we will consider (ordinary) discrete stochastic processes that are sampled versions of the generalized ones: 

s[k] = {s,S(--k)}= s(t)\ t=k 

where s(t) is the continuous-time solution of Q. It should be clear now that the statistics of this discrete process 
are completely specified by Z s (tp) in ( [12) . For instance, we may obtain its ifth-order characteristic function 
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<g{e ji - s ^} = J RK p a (s)e j ^ ds with s = (s[k\, s[k - 1], . . . , s[k - K + 1]) and u = . . . ,Wjr) for any finite 
if by substituting ip = ui±S(-) + W2^(- — 1) + ■ • ■ + wr-<5(- — if + 1) in the characteristic form. Likewise, one can 
determine its correlation sequence by sampling the continuous-time correlation function (as given by ([T] Property 
4]) on the integer grid: R s [kx,k 2 ] = R s (h, t 2 )\ tl=klM=k2 - 

Our objective is now to relate these quantities to the Hermitian-symmetric Green function of the operator LL*. 
The latter, which is the distributional solution of LL*/*^. = can formally be specified as 

where L(oj) (resp., \L(— lS)\ 2 ) is the transfer function of the LSI whitening operator L (resp., LL*). Note that in 
the singular case, the above integral has to be interpreted as a finite part (F.P.) integral in the sense of Hadamard. In 
the event where L _1 is LSI BIBO-stable with impulse response p^, then p-^ L , ~p^*p^. However, in the unstable 
case, the latter convolution product is generally undefined; e.g., poD'(t) = J 71 |]l3p"} (*) = — f 1*1 7^ ( u * uV )W 
where the right-hand side expression is not converging anywhere. Next, we make the link with exponential splines 
by expressing the Green function as a weighted sum of augmented B-splines: 

PLL*(*) = 5>«M/% L .(i-fc) (18) 

fcez 

where q a [k] is the Hermitian-symmetric sequence whose discrete-time Fourier transform is 

Qa(z) = 



\D a (e-^)\i |A a (-o;)|' 
The augmented B-spline kernel /% L , is given by 



'LL* 



fa*$t = A a A> E L. (19) 



where /?l is defined by ( fl5j ). Establishing ( 18 1 is a simple matter of factorization in the Fourier domain. What is not 
so obvious at first sight is that the above entities are always well-defined, irrespective of any stability considerations. 
The generalized exponential B-spline f3-^ L ,(t), in particular, is compactly-supported in [— N, +N] and guaranteed 



to yield a stable expansion (Riesz basis property) |26 Theorem 1]. pjj L . (i) and q a [k], on the other hand, are both 
infinitely-supported; they are either exponentially-decaying (stable scenario with Re(a„) ^ 0) or, at worst, of slow 
(polynomial) growth when n > 0. Our final theoretical tool is a corresponding exponential spline interpolation 
mechanism. 

Property 1 (Exponential spline interpolation): Let f(t) be a function (at most of slow growth) that is included in 



the exponential spline space V^ L » = span{/?L L * (* — k)}k^z C S 1 where /3j^ L , is specified by (19 1 and compactly- 
supported in [— N, N]. Then, 

/(*) = X)/(*)Vi n t(*-*) 

feez 

where ipi nt (x) g Vj; L , is an exponentially-decaying interpolation function whose Fourier-domain expression is 

. /% L » 
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with 

N 

B L (z) = ]T fc h ,(k)z- k . (20) 

k=-N 

Proof: The statement f(t) £ Vj; L , is equivalent to f(t) = Efcez c M/%l* — ^) where C W\ 15 a sequence 
of (possibly slowly-growing) B-spline coefficients. By sampling this expression at the integers and taking the 
z-transform, we obtain F(z) = J2kez f( k ) z ~ k = C(z)B h {z) so that C(z) = F{z)/B L (z). The time-domain 
interpretation is that c[k] = (/i; nt * f)[k] where h lnt is the impulse response of the (inverse) digital filter whose 
frequency response is H- m t(e? u ) — 1/ 'B^e*"). Whenever the purely-imaginary poles of 1/L(u>) are such that 
juj n — jui m 7^ j2wk for any n ^ m and k ^ 0, then /3l generates a Riesz basis J26| Theorem 1], which is 
equivalent to < A < B\ J (eJ u ') < B for any w e M (i and B are the lower and upper Riesz bounds of 
the B-spline basis). Therefore, by Wiener's lemma, we have the guarantee that the sequence h lnt is well-defined 
(/lint € ^i) and exponentially-decreasing because /8tt . (k) is compactly-supported. This leads to the conclusion 
that /(*) = EftezOtat * /)[fc]/% L »(*- fc) = Efcez/[ fc ]^mt(t- fc) where yj int (i) = Efc£Z^mtW/% L ,(i - fc) is 
exponentially-decaying as well. The Fourier transform of this last expression is (p- m t(uj) = Hi nt (e : ' u ')/3^ L ,(uj). ■ 

We are now ready to uncover the relation between the second-order statistical characterizations of the continuous- 
time and discrete-time versions of our stochastic processes. We start with the stationary case where the underlying 
TVth-order system is stable (cf. (Tl Property 1]). 

Property 2 (Conversion from discrete to continuous): Let s be a generalized (Gaussian or non-Gaussian) station- 
ary process that satisfies the iVth-order stochastic differential equation Q with a white noise excitation. Then, the 
correlation functions of the continuous-time and discrete-time (e.g., sampled) instances of the process are linked 
through the interpolation formula 

r s (t) =£{s{t') -s(t' + t)} = £V s [% int (t-fc) 

fcez 

where ipi nt (x) is specified in Property [T] and r s [k] = r s (t)\ t _ k . The Fourier-domain counterpart of this expression 
provides the exact link between the continuous and discrete-domain power spectra of the process: 

Remark on notation: While we are using a common symbol to denote the continuous and discrete autocorrelation 
(resp., power spectrum) of s, we are relying on the index variables to distinguish between the two settings. 
Specifically, $ s (w) = J 7 {r s (t)}(uj) is the Fourier transform of the continuous-time autocorrelation function r s (t), 
while & s (z) is the z-transform of the discrete-time correlation sequence r s [k] (or, equivalently, the discrete-time 
Fourier transform if we set z = e^). 

Proof: Since the discrete process is the sampled version of the continuous one, we have that r s [k] = r s (t)\ t _ k , 

2 

or equivalently, $ s (e JW ) = Enez'M^ + 2?m )- We also know that r s (t) = 0oPll* (*) and ^s( w ) = ^ ^ , as 
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a direct consequence of the innovation model. Putting these elements together, we find that 

\L(-u + 27m)| 2 



where L(u>) is the frequency response of the whitening filter specified by the reciprocal of (|6j. We then use the 
B-spline connection to show the above ratio is well-defined and equal to <Pi n t(to). To that end, we consider the 
Fourier-domain version of ( fT9| ) 

LL |i(-w)| 2 

together with its periodized counterpart Engx ( w + 27m) = ^ — ^^n)^ = ^ L ( e,JW ) (by Poisson's 
summation formula and the 27r -periodicity of D a (e^)). It now suffices to express the right-hand side of (21 1 as 



the ratio of these two entities, which yields the desired result. The main point of this manipulation is that B^e^) 
is guaranteed to be non-vanishing (due to suitable pole-zero cancellations), while it is not necessarily so for the 
denominator of ( |2T| i. ■ 
In principle, it is possible to obtain similar results for the non-stationary processes covered by [1 Property 2], 
but the correlation formulas are more involved and somewhat harder to get to. We will illustrate the concept by 
considering the case of an TVth-order process with a single order of non-stationarity (no = 1 and ax = ju>o). 
Without loss of generality, we will also assume that the boundary condition is enforced at the origin: s(io) = 
with to = 0. 

Property 3 (Correlation structure of processes with one order of non-stationarity): Let s(t) be an .ZVth-order 
generalized Levy process with whitening operator L that meets the boundary condition s(0) =0 imposed by the 
presence of the pole jujQ on the imaginary axis. The non-stationary correlation of the continuous and discrete-domain 
versions of the process is fully specified by 

R.(t, = £{a(t) ■ W)} = Vs(t' - t) - e jUot v s (t') - e-^'vsi-t) 
R a [k, ft'] = £{s[k] ■ s\p]} = v s [k' - ft] - e j " ok v s [k'] - e - ju ° k ' v s [-k] 

where v s (t) = cr 2 PL L . (t) and where v s [k] — v s (t)\ t= k- These latter entities are linked through the exponential 
spline interpolation formula 

v s{t) = Y V s [ft] Vint (* - ft) , 

fcez 

where <pi n t(t) is specified in Property [T] Moreover, v s (t) — 0(\t\) is a function of slow growth whose generalized 
Fourier transform is given by 

o.M = gg ' = v s (en^M 

\L(~uj)\ 2 

where L(ui) is the frequency response of L which exhibits a zero at ajy = j^o- 

Proof: From JT] Property 4], we have that R s (ti,t 2 ) = o- 2 L _1 L _1 *{^(- — ti)}(t 2 ). Since the system has a 
single singularity on the imaginary axis at juio, (lOi implies that L _1 * = TlsiI^ s where Tlsi is a BIBO-stable 
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LSI system whose transfer function is T(uj) = . Using the Fourier-domain formula (35 1 of I* o Vo with 



(pa = 5 and (£q = 1 and the fact that Tlsi is a standard convolution operator, we find that the Fourier transform of 
TlsiI^W- - h)}(t) is given by 

T(w) 



-j(uj + uj ) 

Likewise, we have that L _1 = I W() .5T LSI which, by considering the complex-conjugate counterpart of the formula 

Q for I« ,5. yields 

r p-j^t! _ pjui ti 

L-^L- 1 *{S(--t 1 )}(t)= / - , |T( W )| : 

Jr -j(w + w ) 



n- r — 12 " l r MI 2 ^— (22) 

PEL- (* - tl) ~ e^p^» (t) e-^p LL , (-ti). 



The critical step in this derivation is the evaluation of integral ( |22| ) which, contrary to appearances, is actually non- 
singular, thanks to the presence of the fourth term in the numerator. To make this explicit, we recall that the proper 
specification of the inverse Fourier transform of l/\L(— uj)\ 2 , which has a second-order singularity at cj = — ujq, 
involves a finite-part integral that can be resolved as follows 



2 



= F.P./|TM| 2 ^%^ 



\uj + L0q\ 2 2t: 

where the numerator is corrected by subtracting the first two terms of the Taylor series of e Jwt around uj = —ljq. 



This means that in order to neutralize the singularity in the denominator of (22i, we need to correct the three 
leading exponentials in the numerator by subtracting their first-order development around cj = — cjrj- This is possible 
by rewriting the numerator as follows 

Numerator = e ^( f -*i) - g-^oCt-d) ^ + j( u + uj )(t - h)] 

- e 3Uotl (e JUJt - e~ juot [1 + j(u + u) )t}) 

- e"^ * {e-^ 1 - e^ otl [1 - j(u + uo)h]) 

= e :M f -*i) _ gi^otigj^t _ g-jwotg-i^ti _|_ e -juj a (t-ti) 

which shows that the terms that are required to make the inverse Fourier integral non-singular precisely add up to 
e -.M)(*-ti)_ g 

C. Discrete increment process 

The important point that has been brought out by the above analyses is that the present class of discrete (or 
continuous-time) processes exhibit long-range dependencies due to the infinite support of their autocorrelation 
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function. This behavior is further exacerbated in the non-stationary case where the (asymptotic) decay is linear at 
best. Fortunately, we have seen that there is a simple way to obtain a much better conditioned signal by applying 
the localization operator A a (cf. JT| Property 3]). The good news is that this concept is directly transposable to the 
discrete domain as well, and that it substantially simplifies the statistical characterization of such signals, irrespective 
of any stability considerations. 

Specifically, the discrete generalized increment process of s[k] is defined as: 

N 

s d [n] = A a s(t)| t=n = ^2 d a [m]s[n - m] (23) 

m=0 

where s(t) is a generalized Ath-order stochastic process with whitening operator L and pole vector a = (ai, . . . , ckjv); 
the discrete AR-type filtering coefficients on the right hand side of ( |23] l are given by 

N N 

D a (z) = d a [m]z- m = -e^z- 1 ). (24) 

m— n— 1 



Property 3 (Characterization of discrete increment process): Let Sd[k] be the discrete increment process associ- 
ated with a (possibly non-stationary) Ath-order generalized process whose characteristic functional Z s (<p) is given 



by (12 1 where L 1 * is the adjoint of L 1 specified by (lOi (see also JT] Eq. (20)]). Then, Sd[k] is stationary with 
an Ath-order of dependency: p Sd (s d [k] \{s d [k - m]} me2 +) = p Sd (s d [k] \s d [k - 1], ... , s d [k - (A - 1)] ). The 
characteristic function of its A'th-order joint probability density function p Sd [s d [k], s d [k ~ 1], . . . , s d [k — (K — 1)]) 
is given by 



A" 



p Sd (u 1 ,...,u> K )=Z w 5^w fc $(t-fc+l) (25) 



where /?l is the generalized B-spline defined by (15 I. The autocorrelation sequence of the process is compactly- 
supported: 

r d [k] = S {s d [k'} ■ s d [k' + k]} = a%fc L , (fc) 
where /?l L , (t) = (fih * while its power spectrum is simply 

* ai (e*") = <r 2 B L (en 



where B L (e^) is defined by ([20 1. 

Proof: The result is a consequence of |T1 Properties 3 and 5]. The pointwise specification (characteristic 
function of order K) is obtained by making the substitution cp = cji S (• ) + • • • + ujk 6(- — K +1) in the characteristic 
form Z Sd ((p) — Z w ((3^ * ip). The independence between s d [k] and s d [k'] for any k' such that \k ~ k'\ > N then 
follows from the fact that the corresponding B-splines are non-overlapping (since the support of /3l is of size N). 
Indeed, due to the independence requirement, the general noise functional has the property that Z w (pi + p^) = 
Z w (pi) ■ Z w (p 2 ) whenever <pi and ipi have non-overlapping support, which is synonymous with independence. 
As for the autocorrelation sequence, it is simply the sampled version of the one given in (T] Property 4]. Likewise, 
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the power spectrum, whose generic form is 

= £rg$^|&,(w + 27rn)| 



reduces to the finite sum CTq X)fcL-JV {k)e"^ k , thanks to the compact support of /?f; L » (*)■ H 
We would like to emphasize that the statistical characterization of the discrete increment process in Property 
[3] is complete and that it covers the full class of Gaussian and non-Gaussian stochastic processes specified by 
the generic stochastic differential equation (ffll, including the unstable scenarios which are outside the classical 
theory of stationary processes. A remarkable feature is the omni-presence of the exponential B-spline kernel /3l, 
which has a fundamental role in all aspects of the characterization. For instance, we observe that the argument 



= XfcLi W n/?L ^ + 1) in the noise functional (25 1 actually corresponds to the generic form of a cardinal 



exponential spline with the Fourier variables taking over the role of the B-spline coefficients. Likewise, the correlation 
structure is entirely specified by the integer samples of /3^ L , (t) (the autocorrelation of while the power spectrum 
is proportional to i?L(e 3 "), the so-called discrete B-spline filter, which also enters the definition of the spline 
interpolator in Property [T] 

The link of course is not coincidental. In spline theory, the construction of B-splines is motivated by the desire 
to find the shortest possible basis functions to represent a certain family of spline functions. Here, the introduction 
of the generalized increment process is aimed at producing a derived signal with the simplest possible statistical 
structure; in particular, the shortest dependency distance. The proposed solution is optimal in the sense that it 
achieves the shortest possible order of dependency, as a consequence of the minimal support property of the B- 
spline. The localization sequence d a [k] is obviously not arbitrary; the guiding principle is that A a must have the 
same null space as L such as to annihilate all the long-ranging exponential/polynomial modes of L _1 . Concretely, 
this is achieved by mapping the continuous -domain poles of the system into the discrete-domain zeros of D a (z) 



via the exponential map z — e s (cf. Eq. (24 1); this also implies that the minimal length of d a [k] is N + 1, which 



puts a lower bound of N on the size of the B-spline. 



D. Discrete innovation models 



Given the fact that the discrete processes s[k] and Sd[k] are linked through the difference equation (23 1, it is 
tempting to investigate whether or not it is possible to go one step further and to specify s[k) through a discrete 
ARMA-type model. Ideally, we would like to come up with an equivalent discrete-domain innovation model that 
is easier to exploit numerically than the defining stochastic differential equation To that end, we perform the 
spectral factorization of the discrete B-spline kernel 

N 

B L (z)= P L L*(k)z- k = B+(z)B-(z) (26) 

k=-N 
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where B^(z) = ^2%=o M z fe = ^l( z1 ) specifies a causal finite impulse response (FIR) filter of size N. The 
crucial point for the argument below is that (e^) (or, equivalently B\\(e^ u ) as in Property [lji is non- vanishing, 
which is equivalent to the requirement that /3l generates a valid Riesz basis [26|. 

Property 4 (Stochastic difference equation): The sampled process of order N with parameters (L, a) satisfies 
the discrete ARMA-type whitening equation 

N JV-1 

d a [k]s[k -n]=Y^ b H k }4 k - m ] 



where d a and are defined by (|24j) and (26i, respectively. The driving term e[k] is a discrete stationary white 
noise (white meaning fully decorrelated or with a flat power spectrum). However, e[k] is a valid innovation sequence 
with independent, identically-distributed samples only if the corresponding continuous-domain process is Gaussian, 
or, in full generality (i.e., non-Gaussian case), if it is a first-order Markov or Levy-type process with N = 1. 

Proof: Since \B£ (e ja, )| = \J B^(e^) is non-vanishing and a trigonometric polynomial of whose roots are 
inside the unit circle, we have the guarantee that the inverse filter whose frequency response is fl+ * is causal- 
stable. It follows that $ e (e JW ) = <7q ^^g^")™^ = °o> which proves the first part of the statement. As for the 
second part, we recall that decorrelation is equivalent to independence in the Gaussian case only. In the non-Gaussian 
case, the only way to ensure independence is by restricting ourselves to a first-order process, which results into an 
AR(l)-type equation with e[n] = Sd[n). Indeed, Property[3]implies that, for N = 1, p Sd (sd[k] \{sd[k — m]} me z+) = 
p Sd (s d [k]). This is equivalent to s[k] having the Markov property since p s (s[k] \{s[k — m]} meZ +) = p Sd {sd[k]) 

= Ps (s[k]\s[k-l}). m 

We remark that the correspondence between continuous-time and discrete-time ARMA models is a classical 
result in the theory of Gaussian stationary processes |10|. The present contribution to the topic is: 1) to make the 
connection completely explicit thanks to the introduction of the localization filter D a (z) and the discrete B-spline 
kernel B-^(z), and 2) the extension of the result for the non-stationary and/or non-Gaussian scenarios. 

IV. Numerical generation of stochastic processes 

A. Determination of B-spline s 

The generalized exponential B-splines were introduced in [26] in order to establish a formal link between the 
continuous-time and discrete-time theories of linear systems. These functions are slightly more general than the 
classical ones specified by ( fl"3] >, which are missing "zeros". Since a differential LSI system is characterized by its 
poles a. = (ax, ■ ■ ■ , «jv) an d zeros 7 = (71, . . . , 7^) with M < N, the idea is to associate it with an identifying 
exponential B-spline function: 

{ ( M \ N 1 _ P a n -jw 1 

p ian) ® = t-' n - n ?CJ _ Q (*)• ( 27 > 

l \m=l / n=l J n ) 

Such B-splines can be computed explicitly on a case-by-case basis using the mathematical software described in 



1 26 Appendix A]. The connection with Eq. (15i is = 6m P(a-,f)(t) where the a n and j m are the roots to 
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the polynomial P N (() = ( N + Un-^' 1 + ■ • • + aiC + a and Q M (0 = b M ( M + &m-iC W-1 + --- + b 1 ( + b , 
respectively. We recall that the generalized exponential B-splines have the following remarkable properties: 

• They are causal and compactly-supported in [0, iV]. 

• They are piecewise-exponential polynomials with joining points at the integers. 

• They are Holder continuous of order N — 1 — M, meaning that the function and its derivatives up to order 
N — 1 — M are bounded. 

• They generate valid Riesz bases (by integer shifting) of the spaces of cardinal L-splines iff. a n — a m ^ 
j2irk, k £ 1 for all distinct, purely imaginary poles. 

• They perfectly reproduce the exponential polynomials that are in the null space of the operator L, as well as 
any of its Green functions pi,, which all happen to be special types of L-splines (with a minimum number of 
singularities). 

The basic operations of the corresponding B-spline calculus are: 

• Convolution by concatenation of parameter vectors: (|8{ 0l?r ) * /3^ a2 -^ 2 ))(t) — /3( ai ; a2 ; lin2 )(t) 
. Mirroring by sign change: (ayy )(-t) = (-1) M e Q ") + N) 

• Complex-conjugation: /3/ arf )(t) = /3( a .^)(£) 

• Modulation by parameter shifting: ei Uot f3/ a .^ (t) = ^(a+juo^+juo) (*) w ^ trl trie convention that j = (j, . . . , j). 
It follows that the autocorrelation B-spline /?l L , = /?l * /?l tnat i s cent ral to our formulation is given by 

= b 2 M {-l) M (j] e Q " ] (3 {1 ,:- a ^.-y)(t + N) (28) 

B. Discrete inverse operators 

We have seen that the discrete increment process has a much simpler statistical structure than the process from 
which it is derived. This is not only advantageous for the analysis of such stochastic processes, but also exploitable 




for synthesis purposes. The latter calls for a discrete operator mechanism for inverting the difference equation (23 1. 
The technique that we propose is in all points analogous to the continuous-domain method presented in Section 



II-B The principle is to factorize A a = A ajv • • ■ A Ql where each individual operator A ari actually corresponds to 
a discrete FIR filter with transfer function D an (z) = 1 — e a ™;z -1 . 

Formally, the inverse operator of A Q is the digital filter whose impulse response h a [k] is the inverse z-transform 
of 1 _ J. i = ^ - a z . Classical system theory tells us that such a first-order filter is causal-stable iff. its z-domain 
pole z p — e a is inside the unit circle, which is equivalent to Rc(a) < 0. It is also possible to change the domain of 
stability to Re(a) > by switching to an anti-causal response instead of a causal one. The corresponding definition 
of the impulse response is 

u[k]e ak if Re(a) < 

— u[— k — l]e ak else 



h a [k] = 
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which is the sampled version of p a (t) in JT] Eq. (22)] (if one excludes the point of discontinuity of u(t) at t = 0). 
The critical configuration is Re(a) = in which case h a [k] is still bounded — but not in t\ — meaning that the filter 
is no longer stable. 

At any rate, the main point is that A~*x[fc] = (h an * x)[k], and that these first-order inverse filters can be 
implemented recursively as: 

Causal recursion for Re(a„) < 

y[k] = (h an * x)[k] = e a "y[k - 1] + x[k] 

Anti-causal recursion for Re(a n ) > 

y[k] = (h an * x)[k] = e~ Q " (y[k + 1] - x[k + 1]) 

The final ingredient is the discrete counterpart of the operator l Wa! s specified by (|9j; that is, the unique right inverse 
of A JUo that sets the output signal to zero at k = 0. This operator, which is denoted by A7j o s , is given by 

= (h juo * KM - e j " ok (h juo * x)[0] 

where the second term is a properly- weighted complex sinusoid that is in the null space of A 3 - Wo . For k > ko, the 
above formula simplifies to 



,ju) (k-rn) 



which is an expression that can also be updated recursively. If k < 0, the summation bounds are simply interchanged. 



Using the same notation and pole ordering as in Section II-B we are then able to specify a global right inverse of 
A„ as 



A- 1 



A 



A" 1 ■••A~ 1 , 



(29) 



shift-variant 



LSI part 



which are used to specify the corresponding continuous-domain boundary conditions (111. Let s[k] — A a 1 r[fc] 
where r[k) be an arbitrary input signal. Then, the above operator imposes the tiq boundary conditions 

s[0] = 

WM " (30) 

^ Ajn-'-Aj^WM = 0, 
while its right-inverse property ensures that A a s[/c] = A a A cx 1 r[fc] = r[k\. In the stationary case where uq = 
(i.e., Re(a„) ^ 0, n = 1, . . . , N), we also have that A" 1 A a r[fc] = r[k] (left-inverse property). 

A small word of caution is in order here. The above discrete-domain boundary conditions are only rigorously 



equivalent to the continuous -domain ones in (111 for no < 1. Indeed, it is illusory to attempt imposing exact 
constraints on the derivatives of such signals if all we have at our disposal are samples on a discrete grid. The good 
news, however, is that A JLJn is, by construction, the best first-order approximation of the continuous-domain operator 
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D— ju) n I with the property that : Aj U s(t) — ifij UJn * (D — jui n Id)s) (t) where /3j u is the corresponding first-order 
B-spline. In particular, the latter equation ensures convergence to the exact derivatives as the reconstruction grid 
gets finer (in the same way as finite differences tend to derivatives as the step size goes to zero). 

The theoretical alternative is to accept the discrete-domain boundary conditions as they are, assuming that we can 
properly map them back into the continuous domain. This is indeed feasible by extending our notion of continuous- 
domain boundary conditions, as shown in the appendix. The main point is that there is a unique right inverse of L 
that is admissible (in the sense of (T] Theorem 3]) and compatible with the "discrete" boundary conditions ( [30| : it 
is described in the last paragraph of the appendix. 

C. Algorithms 

1) Gaussian case: The generation of the samples of a generalized Gaussian random process is straightforward 
since we can rely on the equivalent discrete innovation (ARMA) model in Property [4] Given a set of parameters 
a (poles), 7 (zeros), and Cq (noise variance), the procedure is then as follows: 



• Computation of Bi,(z) and spectral factorization as in (26 1. 

• Generation of the innovation signal e[k] which is a random sequence of i.i.d. Gaussian random variables with 
zero mean and variance <Tq. 

• FIR filtering with ar, d inversion of the model via the application of the inverse operator A" 1 which may 
be time-invariant or not, depending on the type of process. 

2) Poisson case: This case is slightly more difficult, but can still be handled exactly by starting from the 
generalized increment process Sd[k], Here, we are using the fact that a realization of a Poisson noise with parameter 
(X;p a (a)) has the explicit form 

w(t) = ^2a n S(t - t n ) 

n 

where t n are random, uniformly-distributed locations over the real line (point process) with an average density of 
A, and where the amplitudes a n are i.i.d. random variables with PDF p a (a)- If we now restrict the observation of 
the process over a time interval [0, T], the generation may proceed as follows: 



Analytical computation of the B-spline /3jj(t) = fejvf P(a;i){t) using formula (27i. 

Generation of the point process (t n ) over the slightly enlarged interval [-N, T+N] together with the amplitude 
variables a n . 

Exact computation of the corresponding discrete increment process by appropriate resampling of the B-spline 
functions: 

s d [k] = (/3 L * w)(t)\ t=k = ^2 a n/3h(k - t n ) 



• Inversion of the model via the application of the inverse operator A^ 1 which, again, may be time-invariant or 
not. 

In effect, the continuous-time realization of the stochastic process s is a non-uniform L-spline with knots at the t n . 
Its explicit analytical form is sit) — poit) + 2~2 n ctnph(t — t n ) where po(t) is a component that is in the null space 
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of L and Ph{t) is a Green function of L. In the stationary scenario, po(t) may be seen as a random component that 
condenses all impulsive noise contributions from outside the generation interval. In the non-stationary case, it has 
the stricter role of enforcing the no boundary conditions imposed by the presence of poles on the imaginary axis. 

3) Alpha-stable case: Here, we can benefit from the key property that any filtered version of an alpha-stable 
noise remains alpha-stable. Indeed, the characteristic function of the variable x — (w, tp) where w is an SaS noise 



(cf. specification of f a (cj) in Section II-C i is given by 



exp I -°-r\\^\\t a 
= exp ,HL' 

where \\tp\\1 = J R |v(i)| a dt is a normalization constant that is shift-invariant; that is, \\ip(- — £o)H2 = IMIl ■ 
This implies that x has an alpha-stable distribution, and by extension, that any linear transformation of an alpha- 
stable process s is alpha-stable as well (23), (27). It is therefore a simple matter to generate an alpha-stable Markov 
process (N — 1) whose increments are independent (cf. Property|4]). The situation gets more delicate for higher-order 
processes because of the necessity of generating an alpha-stable discrete increment sequence with an TVth-order of 
dependency. The first approach that comes to mind is to run an adapted version of the Gaussian algorithm where 
the discrete input noise is alpha-stable instead of Gaussian. Since alpha-stable laws are preserved through linear 
combinations, this will at least ensure that the marginals are alpha-stable and that the second-order dependencies 
are the correct ones. This discrete innovation approach, however, is not entirely satisfactory because decorrelation 
is not rigorously equivalent to independence. 

The alternative approach that we propose is to use a piecewise-constant approximation of the B-spline /?l with 
an oversampling factor of m: 

mN , , 

0L,m(t) = V/3 L (fc/m)A) [m(t ) 

* — ' \ TO 

where /3o(mt) = l[o,i)(£) is a rectangular function of size 1/to. The basic results from approximation theory 
ensure that linim^oo /3L,m(t) = pointwise and in all L p -norms with the error decaying like 1/to (since 

piecewise-constant splines have first-order of approximation). Starting from the oversampled version of the first- 
order alpha-stable increment process Sd : i{k/m) = (/3 (to-) * w) (k/m), which is an i.i.d. alpha-stable sequence, we 
are then able to compute the samples of the discrete increment process by applying the following convolution-like 
equation 

mN ( k \ 

{ph,m * uOWUtf = Yl fcW™<) s d,i I m(k' -—)). 

k=0 ^ m ' 

The approximation can be made arbitrary close by increasing the over-sampling factor to. The computational 
overhead is essentially that of generating to times more i.i.d. random variables as in the Gaussian algorithm. The 
remainder of the procedure is the same as in the Poisson case. Note that this algorithm is generic and applicable 
to other types of Levy noises as well. 
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poles = [0], zeros = [] 




(a) B-spline (solid) and autocorrelation (dashed) (b) Gaussian 




(c) Poisson (A = 1/32, Gaussian amplitudes) (d) SorS (or= 1.2) 



Fig. 2. Example 1: Generation of generalized stochastic processes with whitening operator L = D (pole vector a = (0)): (a) B-spline 
functions = rect(t — g) and /3jj L » (i) = tri(t), (b) Brownian motion, (c) Compound Poisson process with A = 1/32 and Gaussian 

amplitude distribution p a (o) = (2ir)~ 1 / 2 e~ a / 2 , c) SaS Levy motion with a = 1.2. 



We conclude this section by indicating that we can also arbitrarily change the sampling step (which had been 
set to T = 1 for simplicity) via a simple rescaling of the poles, zeros and noise variance. The main point of the 
argument is that L^.^Tus) — T M ~ N L( a /T-,-y/T) ( w ) an d that the white noise property is invariant to dilation (up 
to a normalization factor). 

V. Illustrative Examples 

Examples of realizations of Gaussian versus sparse stochastic processes are shown in Figs. 2 to 5. These signals 
were generated using the algorithms described in Section [TV-C| for the three types of driving noises: Gaussian (panel 
b), impulsive Poisson (panel c), and symmetric-alpha-stable (SaS) with a = 1.2 (panel d). 

The relevant operators are: 

• Example 1: L — D (Levy process) 

• Example 2: L = D 2 (second-order extension of Levy process) 

• Example 3: L = (D — aiId)(D — a^Id) and a = (j37r/4, — j37r/4) (generalized Levy process) 

. Example 4: L = (D - a 1 Id)(D - a 2 Id) and a = (-0.05 + jn/2, -0.05 - j'tt/2) (CAR(2) process) 
The corresponding B-splines (/3l and are snown m the upper left panel of each figure. 

The signals that are displayed side-by-side share the same whitening operator, but they differ in their sparsity 
patterns which come in three flavors: none (Gaussian), finite rate of innovation (Poisson), and heavy-tailed statistics 
(SaS). The Gaussian signals are uniformly textured, while the generalized Poisson ones are piecewise-smooth by 
construction. 
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poles = [0, 0], zeros = [] 





(b) Gaussian 



(c) Poisson (A = 1/32, Gaussian amplitudes) 



(d)SorS (cr= 1.2) 



Fig. 3. Example 2: Generation of generalized stochastic processes with whitening operator L = D 2 (pole vector a. = (0,0)): (a) B-spline 
functions /3l(*) = tri(t) and /3^ L , (t) (cubic B-spline), (b) Gaussian process, (c) generalized Poisson process with A = 1/32 and Gaussian 
amplitude distribution, c) generalized SoS process with a = 1.2. 



A. Self-similar processes 

The classical Levy processes (Fig. 2) are obtained by integration of white Levy noise; they go hand-in-hand 
with the B-spline of degree (rect), and its autocorrelation (triangle function) which is a B-spline de degree 1. 
The Gaussian version (Fig. 2b) is a Brownian motion. It is quite rough and nowhere differentiable in the classical 
sense. Yet, it is mean-square continuous due to the presence of the single pole at the origin. The Poisson version 
(compound Poisson process) is piecewise-constant, each jump corresponding to the occurrence of a Dirac impulse. 
The SaS Levy motion exhibits local fluctuations punctuated by large (but rare) jumps, as is characteristic for this 



type of process |23|, |28'|. Overall, it is the jump behavior that dominates making it even sparser than its Poisson 
counterpart. 

The example in Fig. 3 (second-order extension of a Levy process) corresponds to one more level of integration 
which yields smoother signals (i.e., one-time differentiable in the classical sense). The corresponding Poisson process 
is piecewise-linear, while the SaS version looks globally smoother than the Gaussian one, except for a few sharp 
discontinuities in its slope. The basic B-spline here is a triangle, while /3j L * is a cubic B-spline. The signals in Fig. 
2 and 3 are non-stationary; the underlying processes have the remarkable property of being self-similar (fractals) 
due to the scale-invariance of the pure derivative operators. The Gaussian and SaS stable processes are strictly 
self-similar in the sense that the statistics are preserved through rescaling. By contrast, the scaling of the Poisson 
processes necessitates some corresponding adjustment of the rate parameter A fl6). 
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poles = [j37i/4, -j37r/4], zeros = [] 




(a) B-spline (solid) and autocorrelation (dashed) (b) Gaussian 




(c) Poisson (A = 1/32, Gaussian amplitudes) (d) SorS (or= 1.2) 



Fig. 4. Example 3: Generation of generalized stochastic processes with whitening operator L = (D — aiId)(D — «2ld) and a = 
(j3n/4, —jSn/A): (a) B-spline functions /3l and /%; L *, (b) Gaussian process, (c) generalized Poisson process with A = 1/32 and Gaussian 
amplitude distribution, c) Generalized SaS process with a = 1.2. 



poles = [-.05 + j7i/2, -.05 - jTT/2], zeros = [] 




(a) B-spline (solid) and autocorrelation (dashed) (b) Gaussian 




(c) Poisson (A = 1/32, Gaussian amplitudes) (d) SaS (or = 1.2) 



Fig. 5. Example 4: Generation of generalized stochastic processes with whitening operator L = (D — aiId)(D — cV2ld) and a = 
(—0.05 + jir/2, —0.05 — jTr/2): (a) B-spline functions /3l and /?l L *, (b) Gaussian AR(2) process, (c) Generalized Poisson process with 
A = 1/32 and Gaussian amplitude distribution, c) SaS AR(2) process with a = 1.2. 
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B. Bandpass processes 

The second-order signals in Fig. 4 are are non-stationary as well, but no longer self-similar. They are real-valued, 
and C 1 -continuous almost everywhere (pair of complex-conjugate poles in the left complex plane). They constitute 
some kind of modulated (or bandpass) counterpart of the Levy processes which appears to be much better suited for 
the modeling of acoustic signals. As in the other examples, the Gaussian version is looking cluttered. The Poisson 
signal is somewhat stereotyped (stretches of pure oscillating regime) and not quite as realistic looking as its SaS 
counterpart. 

As soon as the poles are moved away from the imaginary axis, the processes become stationary. This is illustrated 
in Fig. 5 with some CAR(2) (continuous autoregressive) examples, the non-Gaussian versions of which having a 
marked tendency to exhibit characteristic bursts associated with the impulse response of the system. These latter 
processes are part of the stationary CARMA family characterized by Brockwell using an alternative stochastic 
integration/state-space formulation p7| . 

C. Mixed processes 

One can also construct signals with a more complex structure by simple addition of independent elementary 
processes. This results into a mixed process, s m ; x = Si + • • • + sm, whose characteristic form is the product of the 
characteristic forms of the individual constituents: 

M / „ M 

Z Sml M = II = ex P / E /™( L m *¥>(*)) dt 

where s m is some elementary process with whitening operator L m and Levy function f m {uj). As a demonstra- 
tion of concept, we have synthesized some acoustic samples by mixing random signals associated with elemen- 
tary musical notes (pair of poles at the corresponding frequency). These can be downloaded from the web at 
http : / /bigwww. epf 1 . ch/sparse| The Gaussian versions are diffuse, cluttered and boring to listen to. Our 
generalized Poisson and SaS samples are more interesting perceptually — reminiscent of chimes — with the latter 
sounding less dry and more realistic. Note that mixing does not gain us anything in the Gaussian case because the 
resulting signal is still part of the traditional family of Gaussian ARMA processes (this follows from Parseval's 
relation and the fact that J2m=i jX — 1° * s expressible as an equivalent rational power spectrum). This is not 
so for the non-Gaussian members of the family, which are generally not decomposable, meaning that the mixing 
of sparse processes opens up new modeling perspectives. Interestingly, the Gaussian acoustic samples are almost 
impossible to compress using mp3/AAC, while the generalized Poisson and SaS ones can be faithfully reproduced 
at a much lower bit rate. 

VI. Conclusion 

The main point of this paper has been to show that the spline interpretation that links the continuous- and discrete- 
time deterministic linear system theories has a direct counterpart in the linear theory of stochastic processes. While 
the connection between SDEs and stochastic difference equations is well understood in the classical framework of 
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Gaussian stationary processes, it is much less so when (i) the excitation noise is non-Gaussian, and/or (ii) when the 
underlying system is unstable. We have argued that these two extensions are essential for producing signals that 
are sparse — which calls for non-Gaussian excitations — and compressible in a wavelet basis (because self-similar 
processes are solutions of unstable SDEs). Our main effort in this series of papers has been to address these issues by 
setting the foundation of a general framework that extends the bounds of the traditional theory of Gaussian stationary 
processes. The good news is that our generalized formulation leads to a simple universal conversion scheme by 
which a stochastic differential equation is mapped into some corresponding stochastic finite difference equation. 
The cornerstone of this approach is the existence of a compactly supported exponential B-spline, /3l, which acts as 
the mathematical translator between the continuous domain operator L and its discrete version Ld- The elucidation 
of this A-to-D connection has direct implications for signal synthesis (generation of sparse stochastic processes) 
and statistical analysis (proper specification of likelihood functions, optimal signal estimation). Most importantly, 
it provides a functional approach that facilitate the derivation of the joint statistics of such processes, especially in 
the non-Gaussian cases. 

While the proposed framework opens up new modeling perspectives, it also calls for further mathematical 
investigations. In particular, we are currently applying our results to quantify the sparsifying properties of wavelet- 
like expansions where the wavelets are matched to the operator; e.g., ip = L*(f> where <p is an appropriate smoothing 
kernel |29|. We are also postulating that the smoothness properties (Holder and Sobolev exponents) of our extended 
family of CARMA processes are directly related to those of the underlying B-splines. While this is justifyable in 
the Gaussian and Poisson cases fl6) , |30], the details still need to be worked out for the other brands of noise, 
especially the ones with unbounded variance (e.g., SaS) for which a mean-square interpretation cannot be provided. 



The guiding principle for defining non-stationary processes with generalized boundary conditions is to extend 



and where tpo(t) is some given compactly-supported function such that (po(— ojq) 7^ 0. We note that the above 
operator is well-defined pointwise for any / G L\ and that it is a right inverse of (D — jojQld) because the 
sinusoidal correction on the right is in the null space of the operator. By design, Io; , Vo is such that it imposes the 
generalized boundary condition 



Appendix: Generalized boundary conditions 




(31) 




(I. 



f,<Po)=0 



(32) 



for any input function /. 
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Our next task is to show that the adjoint of this operator is admissible. To identify I* o , we perform the 
inner-product manipulation 



(U, Vo f,g) = (Uf,g)~(e^\g) 

= (/.C 5> -g{-u ) 



which, by identification with (/,I£ Vo g), yields 

T^ m m = ^{ g -£=^ )V0 }(t) 03) 

where I* o is the anti-causal convolution operator whose impulse response is pj U (t) = u(— t)e~i Uot . The right- 
inverse property of I WOj¥ , automatically gets transposed into a left-inverse property for its adjoint I* . Next, by 
using the fact that e JU,D * (^f(t) — ;J~~f3^^¥>o(i)J dt = and applying the same technique as in the proof of 
(T] Proposition 2], we show that 

IT* fffll <r r II /lioo,y 

for any / 6 L<x>,r (the space of functions with algebraic decay of order r) where C Vo is a constant that solely 
depends upon ipQ, This proves that I* VQ is a continuous operator on 1Z (the space of rapidly-decreasing functions), 
and, by implication, a continuous map from S into L p with p > 1. The same holds true for any combination 
(iteration) of such elementary operators. 

For completeness, we are giving the equivalent Fourier-based definition of the relevant pair of inverse operators 
which are valid for distributions as well: 

/ A, .+ A, .„+ ihr, ( — \ 

du> 




_ ju t yo(-") 
yo(-"o) 

j(ui - w ) / 2tt 



(34) 



P n/ , / I ' /qVZ'W^M .'art dtJ n rx 

W" = / i -j(a, + W0 ) "J 6 ^tT' (35) 

Observe that both Fourier integrals are non-singular and that we recover the formulas in (T| Table 1], as well as 
(|9j, by setting <^o = £(■ — to) and <ySo = 5, respectively. 

We can now replicate the construction of an admissible left-inverse operator L _1 * for the general TVth-order 
differential system in |T| Section III-B]. In the case of an nrjth-order of singularity, the generic form of a proper 
inverse operator that is admissible in the sense of Q is 

where Tlsi lS some "standard" <S -continuous convolution operator. The adjoint L _1 = I* n l/)n '"^ui ^i^LSI' 
which is the right-inverse of L, is then such that it imposes the generalized boundary conditions on the output 

3 The derivation of the first formula relies on the duality-product version of Parseval's relation: (s, ipo) = L,a(t)tpo(t) dt = 
K Jk*( w )^o(-^) where 0q(-u;) = F{W)}- 
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signal s = L 1 w 



(tp no ,s) = 

(yno-l>( D -J w no Id ) S ) = 

( ( p 1 ,(p-ju 3 Id)--'(p-ju no li)a) = 0, 



(36) 



for any driving term w. 

Interestingly, if we select ip nQ 



n -l 



M (jw„ -i,jw„ )' 



91 



/3/. s, we end up 



with a set of continuous-time boundary conditions (36 1 that is rigorously equivalent to the "discrete" one in (30 1. 



Since the specification of boundary conditions is somewhat arbitrary anyway, this is clearly our preferred choice. 
It has the advantage of ensuring a perfect compatibility between the continuous and discrete-domain specifications 
of these processes. 
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